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Abstract 

A fractional relaxation equation in dielectrics with response function of the Havriliak-Negami 
type is derived. An explicit expression for the fractional operator in this equation is obtained and 
Monte Carlo algorithm for calculation of action of this operator and is constructed. Relaxation 
functions calculated numerically according to this scheme coincide with analytical functions 
obtained earlier by other authors. The algorithm represents a numerical way of calculation in 
relaxation problems with arbitrary initial and boundary conditions. A fractional equation for 
electromagnetic waves in such dielectric media is obtained. Numerical results are in a good 
agreement with experimental data. 
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1. INTRODUCTION 

When a step-wise electric field is applied, polarization of 
a material approaches its equilibrium value not instantly 
but after some time later. Hereditary effect of polarization 
is expressed through the integral relation. For electric in- 
duction in an isotropic medium, we have [Frohlich (1958)] 

t 

D(t) = £oc E(t) + J K(t- t') E(t')dt'. (1) 

— oo 

Fourier transformation of the last expression gives 
D(w) = £*(cj)E(cj), e*(ui) = Eoo + k(ui) 

Relaxation properties of different media (dielectrics, semi- 
conductors, ferromagnetics, and so on) are normally ex- 
pressed in terms of time-domain response function <p(t) 
which represents the current flowing under the action of a 
step-function electric field, or of the frequency-dependent 
real and imaginary components of its Fourier transform: 

oo 

4>{uj) = { e-' lut (j){t)dt = (j>'(ui) - i<j>"{u). 



The classical Debye expression (see Debye (1954)) for a 
system of noninteracting randomly oriented dipoles freely 
floating in a neutral viscous liquid is 

4>{") = t - i 

1 + ILOT 

where r is the temperature-dependent relaxation time 
characterizing the Debye process: 

(j)(t) = 0(O)exp(-i/T-)> t>0. 

The latter function obeys the simple differential equation 

#(t) 



clt 



VW = o. 



Numerous experimental data gathered, for instance, in 
books (Jonscher (1977), Ramakrishnan & Lakshmi 
(1987)) convincingly show that this theory is not able 
to describe relaxation processes in solids, the relaxation 
behavior may deviate considerably from the exponential 
Debye pattern and exhibits a broad distribution of relax- 
ation times. There exist a few other empirical response 
functions for solids: the ColeCole function (Cole & Cole 
(1941)) 



i 



1 + (iu/ujpY 



< a < 1, 



(2) 



The Fourier transform of the response function <f>(t) is 
related to the frequency dependence of dielectric permit- 
tivity through the following relation 



the Cole-Davidson function (Davidson & Cole (1951)) 

^ = n + f 1 / w? ' < /3 < 1, (3) 
(1 + (iuj/uj p ))p 

and the Havriliak-Negami function (Havriliak & Negami 
(1966)) 



Here e s = e* (0) is the stationary dielectric permittivity. 



i 



(i + (iu/uj p )<x)e> 



< a,fi < 1. 



(4) 
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Here u> p is the peak of losses, a and (3 are constant 
parameters. 



In the present paper, we obtain a fractional relaxation 
equation in dielectrics with response function of the 
Havriliak-Ncgami type. With the help of definitions of the 
functional analysis we derive the explicit expression for 
fractional operator in this equation. Then we construct 
a Monte Carlo algorithm for calculation of action of this 
operator and of the inverse one. The algorithm represents 
a numerical way of calculation in relaxation problems 
with arbitrary initial and boundary conditions. Then, we 
consider propagation of electromagnetic waves in such 
dielectric media and discuss memory phenomenon. 

2. H AVRILI AK-NEG AMI' S RELAXATION 

The most general approximation for frequency dependence 
of response function is given by two-parameter formula 
proposed by Havriliak & Negami (1966). The solution of 
the corresponding fractional differential equation 

[1 + (r D t ) a f<Kt) = 5(t), t = 1/w p , (5) 

based on the expansion of fractional power of operator sum 
into infinite Newton's series 



[l + (T D t )Y = E(f) (roD t ) 
n obtained by IN 

_^ oo 



a(P-n) 



has been obtained by Novikov et al. (2005): 

-l)T(n + 0) (t\ a{n+fi) 



The operator [l + ojp a ( D t ) a ] can also be presented in 
the form (Nigmatullin & Ryabov (1997)): 

\Nlff{t) = [l+u-« D«f/(t) = 
w -«*exp (-^ oD, 1 -) „Df cxp oDj-) f(t) 

The HN function is considered as a general expression 
for the universal relaxation law (Jonscher (1996)). This 
universality is observed in dielectric relaxation in dipolar 
and nonpolar materials, conduction in hopping electronic 
semiconductors, conduction in ionic conductors, trapping 
in semiconductors, decay of delayed luminescence, surface 
conduction on insulators, chemical reaction kinetics, me- 
chanical relaxation, magnetic relaxation. Despite of quite 
different intrinsic mechanisms, the processes manifest as- 
tonishing similarity (Jonscher, 1996). The situation seems 
to be similar to diffusion processes. Random movements 
of small pollen grain visible under a microscope, neutrons 
in nuclear reactors, electrons in semiconductors are quite 
different processes from physical point of view but they are 
the same process of Brownian motion from stochastic point 
of view. This analogy stimulates search of an appropriate 
stochastic model for the universal relaxation law. Investi- 
gations of such kind have been carried out in the works 
(Weron & Kotulski (1996), Weron (1991), Nigmatullin 
(1984) Nigmatullin & Ryabov (1997), Glockle (1993), 
Jurlewicz & Weron (2000), Coffey ct al. (2002), Dejardin 
(2003), Aydiner (2005)). 

Weron & Jurlewicz (2005) have shown how to modify 
the random-walk scheme underlying the classical Debye 
response in order to derive the empirical Havriliak-Negami 



function. Moreover, they have derived formulas for simula- 
tion of random variables with probability density function, 
Fourier transform of which is the Havriliak-Negami func- 
tion. These relations contain stable random variables. 

Coffey et al. (2002) reformulated the Debye theory of 
dielectric relaxation of an assembly of polar molecules 
using a fractional noninertial FokkerPlanck equation to 
explain anomalous dielectric relaxation. 

Dejardin (2003) considered the fractional approach to the 
orientational motion of polar molecules acted on by an 
external perturbation. The problem is treated in terms of 
noninertial rotational diffusion (configuration space only) 
which leads to solving a fractional Smoluchowski equation. 
This model is in a good agreement with experimental data 
for the third-order nonlinear dielectric relaxation spectra 
of a ferroelectric liquid crystal. 

Here with the help of relations of the functional analysis for 
fractional powers of operators, we derive a new numerical 
algorithm of solution of fractional equations corresponding 
to the Havriliak-Negami response. This algorithm is based 
on Monte Carlo simulation of one-sided stable random 
variables. 

3. FRACTIONAL OPERATOR CORRESPONDING 
TO THE HAVRILIAK-NEGAMI RESPONSE 

Let the equi-continuous semigroup {T t ;t > 0} of Co-class 
be defined on locally convex linear topological B-space F. 
The infinitesimal generating operator A of the semigroup 
{T t } is defined as 

A/ = lim h-\T h -I) f 

fij.0 

with domain 

D{A) = {/ G F; lim h~ l {T h - I) f exists in F}. 

h\.0 

According to S. Bochner and R. S. Phillips, the operators 



T t , a x = T t x = 



_^ ) J t- 1 / Q 3 ( Q )(t- 1 / Q s) T s x ds, t>0, 
o 

x, t = 



where g^ a \t) is a one-sided stable density, constitute an 
equi-continuous semigroup of Co-class. 

The corresponding infinitesimal operators arc connected 
through the following relation (Iosida (1980)): 

A a x = -(-A) a x, Vx e D(A). 

The infinitesimal operator —[1 +_oo D"] generates the 
semigroup 

T h = e- h e- h - D >, 
According to the Bochncr-Philips relation, the semigroup 
generated by the infinitesimal operator [1+_ 00 D"]' 9 , where 
(3 < 1, has the form 

oo 

Ttf = J t- 1/ V a) (i- 1//3 r)TV/dr = 



J dr e- T t-WgWit-Wr) J t^* (r^u) T u f du 
o o 



Considering this integral as averaging over ensembles of 
stable random variables, we obtain the following relation- 
ship 



f t / = /exp(-t 1 /%)/(^-[t 1 /^ 



1A 



S a 



From this semigroup we can obtain the corresponding 
infinitesimal operator [1 +_oo D"]^. 

To find the inverse operator 

[!+_«, D?]-0 0<a,/3<l, 
we use the relation for potential operator 



T s fds. 



Consequently, 



,00 . 
= / jexp(-t^S p ) f[z-{t^Sp) 1,a s}jdt\. 

* ' 
Here and are one-sided stable random variables with 
characteristic exponents < a < 1 and < j3 < 1. 

Introducing exponentially distributed random variable £/, 
we arrive at 



[1 +_ 00 D?]""/ = ( I f(z- S a e /a ) — < 



9 j8 

°/3 



(3(S/U^ 1 f(z-S a u 1 ^)) 



(6) 



We use this formula to find the solution of fractional 
relaxation equation for arbitrary prehistories of charging- 
discharging process. 

4. FRACTIONAL WAVE EQUATION 

Substituting the Havrilyak-Negami dependence of permit- 
tivity on frequency 

e (uj) = £oc + — — g 
into the Fourier transform of the equation (1), we obtain 



exp 



p /- n 1 

a 



l^expl^-t-ooD^jECt) 



a 

Here special forms of fractional operators arise 

w^/W = [i + <*-ooD«F/(t) = 

w -«"exp(-^ -ooDj—) .^Df exp(^ ^D*"' ) /(/> 
The inverse operator has the form 

vc,-"/(*) = [1 + u- a -ocD^m = 



■ uipP exp 



\ a 



I -ooDrj \f exp 



The following asymptotical relationships take place: 



W^/(i) 



-a/3 



Df /(t), i « 1/ Wp 



/(*) 



(7) 



u(t)/u(0) 
1e1 



1e-1 - 



1e-2 



1e-3 
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wsr^/(t) 



.a/3 



f <C l/wp 



(8) 



Maxwell's equations 



47T. 

rot H = — j 

c 

rot E : 



1 <9D 



ldB 

c~dt 



in combination with the material relations 



D = e 00 E+( £s -e 00 )[l+w- i -ooD?]-' 3 E, B = /iH 
lead to the following wave equation 



^ £oo <9 2 E { M(£a-£oo) ri -i nQ!1 -/5 <9 2 E 



s 2 <9i 2 



-V(div E) - V 2 E 



47T/X 9j 



(9) 



At small times we have (£ <C 1/^p) 

^£oo <9 2 E ;«(e s - £oo) a 



dt 2 



-u. 



+V(div E) - V 2 E = 

At large times we have (t ^> 1/wp) 
^e s <9 2 E ^(e a - £oo) 



p -oo^i 



D 2 r af) E- 



c 2 dt 2 



-V(div E) - V 2 E = 



,D 2+Q E+ 



47T/i dj 



The wave equation presented above is concordant with 
equation obtained by Tarasov (2008-2) from Jonscher's 
universal law. 

5. COMPARISON OF NUMERICAL RESULTS WITH 
EXPERIMENTAL DATA 

We compute charging-discharging curves with the help of 
constructed algorithm for a = 0.998. Charging durations 
are 8 = 10, 7.5, 5, 2.5 seconds. The voltage falls according 
to the Debye exponential law for a long time but after 



u (t)/R, A ^ 



V=200 B 




60 min 
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Figure 2. Experimental discharging curves for paper capacitc 
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Figure 3. Comparison of numerical results (2) with experimental 
results (1) for paper capacitor 8 = 300 , t ss 0, 4 c" 1 . Parameter 
a is equal to 0,998. 3 - exponential function; 4 - difference 
between theoretical curve and exponential function. 
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Figure 4. Experimental discharging curve for electrolyte capacitor 
(0 = 3600 ). 

some crossover moment we observe splitting with respect 
to different values of 9 and transition to non-Debye power 
laws. It seems as if the memory returns to the system 
after some interval of time. Such behavior was called in 



(Uchaikin & Uchaikin (2005)) the regeneration memory 
effect. When a = 1 the relaxation follows the Debye law 
and does not depend on 9: the memory is absent. 

Experimental measurements (see, Uchaikin et al. (2008a)) 
were carried out in the following way. At the hrst time 
capacitor was shunted through the resistor R and am- 
meter. Necessary bias voltage was applied by a power 
source. The capacitor was being charged during time 9. 
After charging time 9 capacitor became shunted again. The 
capacitor consists of technical paper impregnated by oil. 
Its stationary capacitance equals 2 • 10~ 6 F. For current 
limitation the resistor 2 • 10 5 Ohm was chosen. Voltage of 
power source is 200 V. 

Thus, experimental data presented in Uchaikin et al. 
(2008a) and given in Fig. 3 confirm the behavior predicted 
by calculation. Relaxation does not depend on the charg- 
ing prehistory in some time domain but after enough large 
time process becomes dependent on charging history. 

In Fig. 4 the comparison of calculated discharging process 
(points) and experimental data (solid line) is given. It was 
used just one parameter adjusted, a (j3 = 1). Theoretical 
and experimental solutions are in a good agreement. Some 
deviation occurs at large times. 

For charging-discharging process the transition from one 
power decaying law (t~ a ) to another (t~ a ~ r ) is theoreti- 
cally predicted in the case of quite large charging times 9. 
This transition is quite difficult for observation in exper- 
iment, because it occurs at large times. Nevertheless this 
change of exponents in power decaying law has been ex- 
perimentally observed for electrolyte capacitor. Charging 
time was equal to an hour. 

6. CONCLUSION 

Let us summarize what we represent in this article. Frac- 
tional relaxation equation (5) and fractional wave equa- 
tion (9) for dielectrics with the response function of the 
Havriliak-Negami type are considered. The explicit ex- 
pression for the fractional operator in these equations is 
obtained. The Monte Carlo algorithm for calculation of ac- 
tion of this operator and of the inverse one is constructed. 
The algorithm is derived from the Bochner-Phillips re- 
lation for a semigroup generated by a fractional power 
of initial infinitesimal operator. The method is based on 
averaging procedure over ensembles of one-sided stable 
variables. 

Relaxation functions calculated numerically according to 
this scheme coincide with analytical functions obtained 
earlier by other authors. The algorithm represents a nu- 
merical way of calculation in relaxation problems with 
arbitrary initial and boundary conditions. 

The case, when a is close (but not equal) to 1 and /3 = 1, 
is considered in more details. Numerical calculations show 
that in this case does almost not depend on the prehistory 
and approximately coincides with up to some value of 
argument, after which the situation changes. Namely, 
the solutions with different prehistories reveal different 
behavior of the inverse power type on the contrary to the 
exponential behavior in the case. Such phenomena named 
memory regeneration phenomena in previous publications 



of the authors are confirmed by experiments with dielectric 
capacitors. 
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